Stable oscillations of a predator-prey probabilistic 
cellular automaton: a mean-field approach. 



Tania Tome and Kelly C de Carvalho 

Institute de Fi'sica, Universidade de Sao Paulo 
Caixa Postal 66318 

05315-970 Sao Paulo, Sao Paulo, Brazil 
E-mail: ttome@if.usp.br 

Abstract. 

We analyze a probabilistic cellular automaton describing the dynamics of 
coexistence of a predator-prey system. The individuals of each species are localized 
over the sites of a lattice and the local stochastic updating rules are inspired on 
the processes of the Lotka-Volterra model. Two levels of mean-field approximations 
are set up. The simple approximation is equivalent to an extended patch model, 
a simple metapopulation model with patches colonized by prey, patches colonized 
by predators and empty patches. This approximation is capable of describing the 
limited available space for species occupancy. The pair approximation is moreover 
able to describe two types of coexistence of prey and predators: one where population 
densities are constant in time and another displaying self-sustained time-oscillations 
of the population densities. The oscillations are associated with limit cycles and arise 
through a Hopf bifurcation. They are stable against changes in the initial conditions 
and, in this sense, they differ from the Lotka-Volterra cycles which depend on initial 
conditions. In this respect, the present model is biologically more realistic than the 
Lotka-Volterra model. 



PACS numbers: 87.23.Cc, 05.65. +b, 05.70.Ln, 02.50.Ga 
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1. Introduction 

The simplest model exhibiting time-oscillations in a two-component system is the model 
proposed independently by Lotka [USE] and by Volterra [I]. In this model the 
individuals of two species are dispersed over an assumed homogeneous space. It is 
implicitly assumed in this approach that any individual can interact with any other one 
with equal intensity implying that their positions are not taken into account. The time 
evolution of the densities of the two species in the Lotka- Volterra model is given by a set 
of two ordinary differential equations [5j El [7J E] and is set up in analogy with the laws 
of mass-action. Depending on the level of description wanted, the approach based on 
mass-action laws, contained on the Lotka- Volterra model, suffices. However, there are 
situations in which the coexistence takes place in a spatially heterogeneous habitat such 
that the population densities can be very low in some regions. In this case we need to 
proceed beyond the mass-law equations and consider the space structure of the habitat. 
In other words, it becomes necessary to analyze the coexistence by taking explicitly into 
account spatial structured models. 

In fact, the role of space in the description of population biology problems has been 
recognized by several authors in the last years P HQl HH H21 H31 HH UHl HE! HZl HH HH1 
[2QI EH [221 [231 [2H ESI [26]. In a very clear manner, Durrett and Levin [11] have pointed 
out that the modelling of population dynamics systems which are spatially distributed 
by interacting particle systems [TTJ [271 12EJ [29] is the appropriate theoretical approach 
that is able to give the more complete description of the problem. We include in this 
approach probabilistic cellular automata (PCA) [291 EH EI] , which will concern us here. 
We refer to interacting particle systems and PCA as stochastic lattice models. They are 
both Markovian processes defined by discrete stochastic variables residing on the sites 
of a lattice; the former being a continuous time process and the latter a discrete time 
process. 

In the present work we study the coexistence and the emergence of stable self- 
sustained oscillations in a predator-prey system by considering a PCA previously studied 
by numerical simulations [25 26J. This PCA is defined by local rules, similar to the 
ones of the contact process [27], that are capable of describing the interaction between 
prey and predator. Here, we focus on the analysis of the PCA by means of dynamic 
mean- field approximations [10], [281 [221 EH1 E21 [33]. In this approach the equations for 
the time evolution of correlations of various orders are truncated at a certain level and 
high order correlations of sites are written in terms of small order correlations. The 
simplest approximation is the one in which all correlations are written in terms of one- 
site correlation, called simple approximation. In a more sophisticated approximation, 
called pair approximation [TUl [29] , any correlation is written in terms of one-site and 
two-site correlations. 

The simple mean-field approximation is capable of predicting coexistence of 
individuals in a stationary state where the densities of each species, and of empty 
sites, are constant. However, it is not capable of predicting possible time oscillating 
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behavior of the population densities and we have proceed to the next order of mean- 
field approximation. The simple approximation, on the other hand, can be placed in an 
explicit correspondence with a patch model [3 Q21 Elj where unoccupied patches can be 
colonized by prey and patches occupied by prey can be colonized by predators that in 
turn may become extinct. In this approximation the PCA can be seen as an extended 
version of the Lotka-Volterra model which includes an extra logistic term related to the 
empty sites. 

The pair-mean field approximation is able to predict possible time oscillating 
behavior of the population densities that are self- sustained and are attained thorough 
Hopf bifurcations. This is in contrast with the Lotka-Volterra model which presents 
no stable oscillations but exhibits instead infinite cycles that are associated to different 
initial conditions. However, from the biological point of view, one does not expect that a 
small variation in the initial densities of prey and predator result in different amplitudes 
of oscillations. Within our approach, a PCA treated in the pair-approximation, the 
oscillations are associated with limit cycles what mean to say that they are stable 
against the changes in the initial conditions. According to our point of view, the pair- 
approximation, in which the correlation between neighboring sites are treated exactly, 
provides a basic description of the predator-prey spatial interactions. For this reason, 
we will refer to the PCA in this approximation as a quasi-spatial-structured model. 

2. Model 

2.1. Probabilistic cellular automaton 

We consider interacting particles living on the sites of a lattice and evolving in time 
according to Markovian local rules. The lattice is the geometrical object that plays the 
role of the spatial region occupied by particles, in a general case, or by individuals of each 
species in the present case. The lattice sites are the possible locations for the individuals. 
Each site can be either empty or occupied by one individual of different species and a 
stochastic variable rji is introduced to describe the state of each site at a given instant 
of time. The state of the entire system is denoted by rj = (r/i, . . . , % . . . , t/jv) where 
A?" is the total number of sites. The transition between the states is governed by the 
interactions between neighbor sites in the lattice and by a synchronous dynamics. 

The probability P^\rj) of configuration 77 at time step I evolves according to the 
Markov chain equation 



where the summation is over all the microscopic configurations of the system, and 
W(r]\i]') is the conditional transition probability from state 77' at time £ to state r\ at 
time i + 1. This transition probability does not depend on time and contains all the 
information about the dynamics of the system. Taking into account that all the sites 
are simultaneously updated, which is the fundamental property of a PCA, the transition 
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Figure 1. Transitions of the predator- prey model. The three states are: prey or 
herbivorous (H or 1), predator (P or 2) and empty or vegetable (V or 0). The allowed 
transitions obey the cyclic order shown. 



probability can be factorized and written in the form [291 EH] 

TV 

W{Ti\rf) = Uv>i{rk\r/), (2) 
i=i 

where Wi(rji\rj') is the conditional transition probability that site i takes the state rji 
given that the whole system is in state rf. Being a probability distribution, the quantity 
w i(Vi\v') m ust satisfy the following properties: Wi(r]i,T]') > and 

5>ifab/) = i- (3) 

Vi 

The average of any state function F(rj) is evaluated by 

(F( V )) e = J2F( V )P {e \v)- (4) 

v 

The time evolution equation for (F(rj)) is obtained from definition (j3J) and equation (pQ). 
For example, we can derive the equations for the time evolution of densities and two-site 
correlations. 



2.2. Predator-prey probabilistic cellular automaton 

To model a predator-prey system by a PCA, the stochastic variable rji associated to site 
i will represent the occupancy of the site by one prey, or the occupancy by one predator 
or the vacancy (a site devoid of any individual). The variable rji is assumed to take 
the value 0, 1, or 2, according to whether the site is empty (V), occupied by a prey 
individual (H) or by a predator (P), respectively. That is, 

f 0, empty (V), 

rn=l 1, prey(H), (5) 
[ 2, predator (P), 

which defines a three-state per site PCA. 

The stochastic rules, embodied in the transition rate Wi{r}i\rf), are set up according 
to the following assumptions, (a) The space is homogeneous, which means to say that 
no region of the space will be privileged against the others, that is, in principle the 
individuals have the same conditions of surveillance in any space region, (b) The space 
is isotropic, which means to say that there is no preferential direction in this space for 
any interaction, (c) The allowed transitions between states are only the ones that obey 
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the cyclic order shown in figured! Prey can only born in empty sites; prey can give place 
to a predator, in a process where a prey individual dies and a predator is instantaneously 
born; finally a predator can die leaving an empty site. The empty sites are places where 
prey can proliferate and can be seen as the resource for prey surveillance. The death of 
predators complete this cycle, reintegrating to the system the resources for prey. 

The predator-prey PCA has three parameters: a, the probability of birth of prey, b, 
the probability of birth of predator and death of prey, and c, the probability of predator 
death. Two of the process are catalytic: the occupancy of a site by prey or by a predator 
is conditioned, respectively, to the existence of prey or predator in the neighborhood of 
the site. The third reaction, where predator dies, is spontaneous, that is, it occurs, with 
probability c, independently of the neighbors of the site. We assume that 

a + b + c=l, (6) 

with < a, b, c < 1. 

The transition probabilities of the predator-prey PCA are described in what follows: 

(a) If a site i is empty, rji = 0, and there is at least one prey in its first neighborhood 
there is a favorable condition for the birth of a new prey. The probability of site % being 
occupied in next time step by a prey is proportional to the parameter a and to the 
number of prey that are in the first neighborhood of the empty site. 

(b) If a site is occupied by a prey, 7^ = 1, and there is at least one predator in its 
first neighborhood then the site has a probability of being occupied by a new predator 
in the next instant of time. In this process the prey dies instantaneously. The transition 
probability is proportional to the parameter b and the number of predators in first 
neighborhood of the site. 

(c) If site i is occupied by a predator, rji = 2, it dies with probability c. 

The transition probabilities associated to the three processes above mentioned can 
be summarized as follows: 

wMv) = cS( Vt , 2) + [1 - fi(v)]S(Vi, 0), (7) 



wi(i\ v ) = fiMSfa, o) + [i - <&(»?)]%, i), (8) 

Wi (2\ V ) = 9i( V )5( Vi , 1) + (1 - c)5( Vl , 2), (9) 

where 

fM = ^£*(%,l), 9i(v) = b jY,5(Vk,2), (10) 

k k 

and the summation is over the four nearest neighbors of site i in a regular square lattice. 
The notation 5(x, y) stands for the Kronecker delta function. These stochastic local 
rules, when inserted in equation (J2J), define the dynamics of the PCA for a predator- 
prey system. 

The present stochastic dynamics predicts the existence of states, called absorbing 
states, in which the system becomes trapped. Once the system has entered such a state 
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it cannot escape from it anymore remaining there forever. There are two absorbing 
states. One of them is the empty lattice. Since the predator death is spontaneous, a 
configuration where just predators are present is not stationary. This situation happens 
whenever the prey have been extinct. In this case the predator cannot reproduce 
anymore and also get extinct, leaving the entire lattice with empty sites. The other 
absorbing state is the lattice full of prey. This situation occurs if there are few predators 
and they become extinct. The remaining prey will then reproduce without predation 
filling up the whole lattice. The existence of absorbing stationary states is an evidence 
of the irreversible character of the model or, in other words, of the lack of detailed 
balance [29J . However, the most interesting states, the ones that we are concerned with 
in the present study, are the active states characterized by the coexistence of prey and 
predators. 



2. 3. Time evolution equations for state functions 

We start by defining the densities, which are the one-site correlations, and the two-site 
correlations. These quantities will be useful in our mean-field analysis to be developed 
below. The density of prey, predator, and empty sites at time step t are defined thought 
the expressions 

pf(l) = (<%,1)>,, (11) 
P t (£) (2) = (*fa,2)>* (12) 

PP(0) = (5( Vi ,Q))i- (13) 

The evolution equations for the above densities are obtained from their definitions as 
state functions, as given by equation (TJ]), and by using the evolution equation for P^(rj), 
given by equation (CQ). The resulting equations can be formally written as 

if +1) (l) = Ml|i7)>* (14) 

pf +1) (2) = ( Wi (2|77)>* (15) 

P/ m) (0) = ( Wi (0\rj)) e , (16) 

where the transition probabilities for this model are given in equations (JTj), (jHJ) and (J9j). 

The correlation between a prey localized at site i and a predator localized at site j 
at time step £ is defined by 

P,?(l,2) = (<%,1)<%,2)>,. (17) 

The other two-site correlations are defined similarly. The time evolution equation for 
the correlation of two neighbor sites i and j, one being occupied by a prey and the other 
by a predator, is given by 

Pg +1 \l,2) H^(%K-(2|77)>*. (18) 

The other two-site evolution equations are given by similar formal expressions. We 
can also derive equations for three-site correlations. Since we are interested here on 
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approximations in which only the one-site and two-site correlations should be treated 
exactly, the above equations suffice. 

We call the attention to the fact that equation (1TB]) includes the product of two 
transition probabilities. This is a consequence of the synchronous update of the PCA 
which allows that both neighboring sites i and j have their states changed at same 
time step. This situation does not occur when we consider a continuous time one- 
site dynamics. Therefore, although local interaction in the present PCA and in the 
continuous time model considered in reference pJJj are the same, the predator-prey 
system evolves according to different global dynamics which leads to different time 
evolution equations for the densities and the correlations. 

The exact evolution equations for the one-site correlations are 



c 



c 



(19) 



W = 7Efji(12) + (l-« 

> i 



(20) 



where the summation in j is over the ( nearest neighbors of site %. To simplify notation 
we are using unprimed and primed quantities to refer to quantities taken at time £ and 
t + 1, respectively. 

The exact evolution equations for the correlations of two nearest neighbor sites j 
and k are 



c 



p, fcn (ooi) - - E ^Wiooi) 



C 



a. 



Pjk( Q1 ) " 7 E Pjkn(012) 



c 



P ijk {101) - - E ^ii*n(1012) 



+ C 



;i - j)p jk (2i) - - E PM 212 ) 



C 



+ 7C E P ifcn(201), 



P jfc n(012) 



E iW 1012 ) 



ij-fc n (112) — 7 X] ^»i*n(2112) 



(1 - 7 )P iA (12) - - E ^(212) 



c 



c 



(21) 



+ 7 (l-c)E^( 1 02), 



(22) 
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and 



p i'k(W) = 7 E 



+ (1 - c) 



(1 - -)P jkn (012) - 7 E ^"(1012) 



C 



c, 



cP jk {22) + P jfc (02) - - £ Pvk(l02) 



c, 



6 

+ r 



P J - Jfc (21) + 2 PM 212 ) > ( 23 ) 

where the summation in i is over the nearest neighbors of j and the summation in n is 
over the nearest neighbors of k. 



3. Mean- field approximation 

3.1. One and two site approximations 

The evolution equation for a density in any interacting particle system which evolves in 
time according to local interaction rules always contains terms related to the correlations 
between neighbor sites in a lattice. The evolution equations for the correlations of two 
neighbor sites includes the correlation of clusters of three or more sites in the lattice and 
so on. In this way we can have an infinite set of coupled equations for the correlations 
which is equivalent to the evolution equation for the probability p( e \rj), described in 
equation ([1]) for the automaton. The scope of the dynamic mean-field approximation 
consists in the truncation of this infinite set of coupled equations [3Q1 [311 [32j [33] . 

The lowest order dynamic mean-field approximation is the one where the probability 
of a given cluster is written as the product of the probabilities of each site. That is, all 
the correlations between sites in the cluster are neglected. For example, let us consider 
the cluster constituted by a center (C) site and its first neighboring sites to the north 
(N), south (S), east (E) and west (W) as shown in figure El 

Within the one-site approximation the probability P(N, E, W, S, C) corresponding 
to the cluster shown in figure [2] is approximated by 

P(N, E, W, S, C) = P(N)P(E)P(W)P(S)P(C), (24) 

where P(X), X = N, E,W, S,C are the one-site probabilities corresponding to each 
site. For some stochastic dynamics models this approximation is able to give qualitative 
results that are in agreement with the expected results. 

In order to get a better approximation we must include fluctuations. The 
simplest mean-field approximation that includes correlations is the pair-mean field 
approximation. This approximation is better explained by taking again, as an example, 
the cluster constituted by a center site which and its four nearest neighbors, shown 
above. Within the pair-approximation the conditional probability P(N, E, W, S | C) is 
approximated by 

P(N, E, W, S | C) = P{N | C)P(E, | C)P(W | C)P{S \ C), (25) 



Stable oscillations of a predator-prey probabilistic cellular automaton 

N 



9 



W 



Figure 2. A site (C) of the square lattice and its four nearest neighbor sites (N, E, 
W, S). 



that is, the conditional probability P(N, E, W, S \ C) is written in terms of the product 
of the conditional probabilities P(X\C), X = N, E,W, S. Now using the definition of 
conditional probability we have 

P(N, E, W, S, C) P{N, C) P{E, C) P{W, C) P(S, C) 



P(C) " P(C) P(C) P(C) P(C) ' 



(26) 



or 



P(JV , E , Wt s , c) = PiWPMWMW (27) 

We see that the resulting probability is written as a function of two-site correlations 
P(X,C), and the one-site correlation P(C). 

3.2. Patch model 

The simple mean-field approximation of the predator-prey PCA describes exactly 
the same properties of an extended Levins patch model [H [M]. That is, the PCA 
with local rules similar to the contact process becomes, in the simple mean-field 
approximation, analogous to the Levins model for metapopulation with empty patches, 
patches colonized by prey and patches colonized by predators. 

In the one-site mean-field approximation we consider that the probability of any 
cluster of sites can be written as the product of the probabilities of each site, as in 
equation (1241) . Using this approach, and writing x = Pi(l), y = -Pi (2), and z = Pi(0) it 
can be seen that the set of equations can be reduced to the following two-dimensional 
map [26] 

x' = x + axz — bxy, (28) 
which is an evolution equation for prey density x, and 

y ' = y + bxy - cy, (29) 
which is an evolution equation for predator density y#. Notice that 

z = 1 — x — y. (30) 

The fixed point of this map are those that represent the stationary solutions x' = x 
and y' = y, and they correspond to the three following solutions x\ = 0, y\ — 0, and 
x 2 = 1, 2/2 = 0, and x 3 = a/b, y 3 = (1 — c/fe)/(l + b/a). The first solution corresponds 
to an absorbing states where both species have been extinct. The second corresponds 
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Figure 3. Phase diagram of the patch model. The continuous line represents 
the transition, c\(p), between the prey absorbing (A) state and the active species 
coexistence (C) state. The dashed line separates the two asymptotic time behavior of 
the active state. 

to an absorbing state where predators have extinct. The third solution corresponds to 

an active state where prey and predator coexist. 

Due to the constraint ([6]), the parameters a, b and c are not all independent and 

only two can be chosen as independent. For this reason it is convenient to introduce the 

following parametrization [10J 

1 — c 1 — c 

a = ^ P> b=^—+p, (31) 

and consider p and c as the independent variables. The parameter p is such that 
— 1/2 < p < 1/2 and < c < 1 as before. This parametrization will useful in the 
determination of the different phases displayed by the model. 

A linear stability analysis reveals that solution the (xi,yi) is a hyperbolic saddle 
point for any set of the parameters a, b and c and so it is always unstable. The empty 
absorbing state will never be reached. A linear stability analysis also shows that the 
solution (x2,y2) is a stable node in the following region of the phase diagram c > c\ 
where 

c 1 (p) = i(l + 2p). (32) 

The active solution is stable in the region c < c\ and is attained in two ways: by an 
asymptotic stable focus, where the successive interactions of the map show damped 
oscillations; or trough an asymptotic stable node. In the phase diagram of figure [3] we 
show the transition line between the prey absorbing state and the active state given by 

c = C\. 

In figure H] it is shown the behavior of the densities against the parameter c, the 
probability of predators death, for the special case p = 0.2. In terms of phase transitions 
what happens is that in the phase diagram there is a transition line separating the 
absorbing prey phase and the active phase which is characterized by constant and 
nonzero densities of prey and predator. 
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Figure 4. Densities of predator and prey as functions of the parameter c for p = 0.2, 
for the patch model. 



We may conclude that the mean-field approximation for the predator-prey 
probabilistic cellular automaton with rules (JTj), dU), and (jHJ) is capable to show, under a 
robust set of control parameters, that prey and predators can coexist without extinction. 
However the map defined by equations fl28|) and (T29]) is not able to describe self-sustained 
oscillations of species population densities. 

3.3. Quasi-spatial model 

In order to find if oscillations in the species populations can be described within a mean- 
field approach we consider a more sophisticated approximation, the pair-approximation, 
where correlations of two neighbor sites are included in the time evolution equations for 
the densities. This is the lowest order mean-field approximation which takes into account 
the spatial localization of neighboring individuals. 

In this analysis we will maintain the correlations of one site and the correlations 
of two-sites in the equations. Correlations of three and four neighbor sites will be 
approximated by means of equation ( 1271) . With these approximation the model is 
described by the following set of five coupled equations 

x' = au — bv + x, (33) 
y' = bv + (l- c)y, (34) 

, ,qu qu 2 , r u w. 

u = aa\ aa — —\ + \{±—pa) — aa—\\u — ab — 

Z Z 2 Z X 

2 

WU V 

+ aac h c[(l — f3b)v — ab — ], (35) 

. irn uv u 2 v, _ ,wu 
v = ab\pa h aa + aall — c) — 

X ZX z 

2 2 

+ ab[— - aiA-] + (1 - c)[(l - (3b)v - ab—], (36) 

XX 2 X 
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Figure 5. Phase diagram of the quasi-spatial model. The upper continuous 
line represents the transition, ci(p), between the prey absorbing (A) state and the 
nonoscillating coexistence (Cno) state. The lower continuous line represents the 
transition, c 2 (p), between the nonoscillating coexistence and the oscillating (Cos) 
coexistence state. The dashed line separates the two asymptotic time behavior of 
the nonoscillating coexistence state. 



and 

/ ,r/, n \ uv u2v ^ /-. \r uw i 

w = ao (1 — pa) aa + (1 — c)[w — aa — 

x zx z 

v 2 

+ c\f3hv + ah— ] +c(l - c)s, (37) 
x 

where a and j3 are numerical fractions defined by a = — 1)/C an d /3 — 1/C where 
( is the coordination number of the lattice. For the present case of a square lattice, 
( = 4 so that a = 3/4 and (3 = 1/4. We are using the following notation: u = P(0, 1), 
v = P(l,2), and w = P(0,2) and also r = P(l,l), q = P(0,0) and s = P(2,2). The 
last three correlations are not independent but are related to others by 

r = x — u — v, (38) 

q = z — u — w, (39) 

s = y — v — w. (40) 

We used the properties P(l, 0) = P(0, 1), P(l, 2) = P(2, 1) and P(2, 0) = P(0, 2), that 
follows from the assumption that space is isotropic and homogeneous. 

We have analyzed numerically the five- dimensional map, described by the set of 
equations ( |33l) . ( |34|) . ( |35l) . ( |36l) and ( |37|) . and we have obtained four types of solutions. 
Two solutions are trivial and are given byx = y = u = v = w = and x = 1, 
y = u = v = w = 0. They correspond to the empty and prey absorbing states, 
respectively. The empty absorbing state, where both species have been extinct is 
an unstable solution and never occurs. However, the prey absorbing state is one of 
the possible stable stationary solutions and is stable above the critical transition line 
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0.1 0.2 0.3 0.4 0.5 0.6 

C 

Figure 6. Densities of predator and prey as functions of c for the quasi-spatial model, 
for p = 0.2. 

c = Ci (p) shown in figure [5j Below this line it becomes unstable giving rise to the active 
state. 

The other solutions correspond to the active states where both prey and predator 
coexist. These solutions are of two kinds: a stationary solution where there is a 
coexistence of the two species with densities constant in time, which we call the 
nonoscillating (NO) active state; and another solution where both population densities 
oscillate in time. This solution corresponds to a self-sustained oscillation of the predator- 
prey system and will be called the oscillating (O) active state. In the phase diagram of 
figure [5] there is a line c = c^ip) that separates the NO and O active phases. Figure [6] 
shows the behavior of the densities as a function of c for p = 0. 

3.4- Oscillatory behavior 

In figure [7] we show an example of self-sustained oscillations of the densities of prey 
and predators as functions of time. The oscillating solutions are attained from the 
nonoscillating solutions by a Hopf bifurcation. The fixed point associated to this solution 
is an unstable center which produces a stable limit cycle as trajectories in the phase- 
space of the predator density versus prey density, as can be seen in figure Notice 
that the oscillations are not damped and have a well defined period which is the same 
for the prey density and for the predator density, which implies that the oscillations are 
coupled. A maximum of predators always follow a maximum of prey. This means that 
the abundance of prey is a condition that favors the increase in the number of predators. 
As the predator number increases the prey population decays. The evanescence of prey 
is followed by a decrease in the predator number, giving conditions for the increase of 
prey population until the cycle starts again. 

A well defined oscillatory behavior is found for many biological population, the 
most famous being the one related to the time oscillations of the population of lynx 
and snowshoe hare in Canada for which data were collected for a long period of time 
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Figure 7. (a) Densities of predator and prey as functions of time and (b) density of 
predator versus density of prey, for the quasi-spatial model, for p = and c = 0.016. 



[3, [H]- If the hare population cycles are mainly governed by the lynx cycle then the 
oscillations shown by the present model reproduces qualitatively some of the features of 
this predator-prey dynamics. 

Next we analyze the behavior of the frequency and amplitude of oscillations. Fixing 
the parameter p and varying the parameter c, we verify that in all the oscillating region 
the frequency of oscillation is proportional to parameter c, 

uj ~ c, (41) 

as can be seen in figure [HJ Low frequencies are associated to low values of c; what means 
that, for small values of c, the greater the predator lifetime the greater will be period of 
the oscillation. As to the amplitude A of the oscillations, we have verified, that fixing 
the value of p and varying the parameter c, it increases as c decreases. Our results show 
that, 

A~( c -c 2 )y\ (42) 

as expected for a Hopf bifurcation and shown in figure El The transition line c = c 2 
from the oscillating phase to the nonoscillating phase can either be obtained by using the 
criterion given by equation (142p or by analyzing the eigenvalues associated to the map 
given by the set of equations (133|) . (134"|) . (1351) . (1361) and (|37|) . This last criterion means 
to find the points of phase diagram such that the real part of the dominant complex 
eigenvalue equals 1. 



4. Discussion and conclusion 



The main result coming from the pair mean-field approximation applied to the predator- 
prey PCA is that it is possible to describe coexistence and self-sustained time oscillations. 
Moreover, these are stable oscillations. Given a set of parameters, just one limit cycle 
is achieved, no matter what the initial conditions are. This property is essential in 
describing a biological system since a small variation in the initial condition can not 
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to 




Figure 8. (a) Frequency of oscillations to versus the parameter c. The frequency 
vanishes linearly as one approaches c = 0. (b) Amplitude A of oscillations versus c 
near the Hopf bifurcation point C2 = 0.019. The quantity A 2 vanishes linearly when 
c — > c 2 in accordance with a Hopf bifurcation. 

modify the amplitude, frequency and mean value of the time oscillation densities of a 
predator-prey system. Similar results were obtained from a continuous time version of 
the present model [10J. Although the simple mean-field equations are essentially the 
same in both versions this is not the case concerning the pair mean-field approximation. 
The time evolutions of the pair correlations for the PCA, presented here, depend on 
higher order correlations (up to fourth) when compared to the ones of the continuous 
version (up to third). 

The model studied here is a spatial structured model with individuals residing in 
sites of a lattice and described by discrete dynamic variables. When we perform simple 
mean-field approximation we neglect all the correlations of sites in the lattice. But we 
take into account that there are limited resources for the surveillance of each species. 
For example in the time evolution equation for the density of prey we have an explicit 
term relative to reaction of birth of prey which is the product of the density of prey x by 
the density of empty sites z = (1 — x — y). This coincides with an extended patch model 
approach for predator-prey systems. The presence of this term is what differs the simple 
mean-field equations from the Lotka-Volterra equations. However, taking into account 
the limitation of space and resources the simple mean-field equations are not sufficient 
to get self-sustained oscillations although able to describe damped time oscillations of 
population densities. 

To get self-sustained time oscillations we had to proceed to the next level of 
approximation in which a pair of nearest neighbor sites is treated exactly. This 
approximation can be seen as representing a pair of nearest neighbor sites immersed 
in a mean field produced by the rest of the lattice. The most important feature being 
the fact that the two sites of this pair can be seen as localized in space. The set of 
five equations which results from the pair approximation for the PCA is indeed able 
to produce self-sustained oscillations of population densities. It presents an important 



Stable oscillations of a predator-prey probabilistic cellular automaton 



16 



property that the Lotka-Volterra model lacks, namely, the oscillating solutions are stable 
and are unique for a given set of the control parameters. 
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